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Abstract 

A  revised  version  of  an  earlier  attempt  to  numerically  solve  Miller's  equations 
for  the  RIGIDICE  model  of  frost  heave  is  presented  that  corrects  earlier  mistakes 
and  incorporates  recent  improvements  in  the  scaling  factors  of  ground  freezing. 
The  new  version  of  the  computer  code  also  follows  the  concepts  of  Object 
Oriented  Numerics  (OON),  which  allow  for  easy  modification  and 
enhancements.  Analysis  of  the  program  is  accomplished  with  the  symbolic 
math  program  MathCad.  A  brief  sensitivity  analysis  of  the  input  variables 
indicates  that  those  parameters  that  calculate  the  hydraulic  conductivity 
have  the  greatest  influence  on  the  variability  of  predicted  heaving  pressure. 


For  conversion  of  SI  metric  units  to  U.S./British  customary  units  of  measurement 
consult ASTM Standard  E380-89a,  Standard PracticeforUseoffheInternational 
System  of  Units,  published  by  the  American  Society  for  Testing  and  Materials, 
1916  Race  St.,  Philadelphia,  Pa.  19103. 
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NOMENCLATURE 


Symbol 

Units 

Definition 

G,W,I 

m“^ 

volume  fractions  of  grains,  water  and  ice 

Wn,  Wn+i 

m~^ 

volumetric  water  content  within  layers  n  and  n+1 
in  the  fringe 

e 

porosity 

Wsao  Wd 

m3  m"^ 

volumetric  water  contents  at  saturation  and  at  the 
lower  limit  of  freezing 

W(,If 

m^m“^ 

volumetric  water  content  and  ice  content  of  the 
frozen  soil  existing  between  two  mature  ice  lenses 

ms”^ 

heave  rate 

m  s“^ 

penetration  rate  of  freezing  front 

flux  of  water 

(flw)b 

m”^  s"^ 

flux  of  water  into  the  fringe 

(^w)n'  (‘7w)n+l 

-2  -1 
m'^m  ^ 

flux  of  water  into  a  layer  in  the  fringe 

m  s“^ 

hydraulic  conductivity 

(^w)sat 

m  s“^ 

saturated  hydraulic  conductivity 

a,  p 

— 

Brooks  and  Corey  coefficients 

Mw/  “i 

Pa 

water  and  ice  gauge  pressures 

<t> 

Pa 

capillary  pressure  {u^-Uy^) 

<t'b 

Pa 

capillary  pressure  at  base  of  fringe  (i.e.,  ice  entry 
pressure) 

V<t>j 

Pa 

capillary  pressures  within  layers  n  and  j  in  the 
fringe 

Pa 

total,  effective  and  neutral  stresses  as  expressed 
by  the  Terzaghi  equation 

X 

— 

Snyder-Bishop  stress  partition  factor 

fw 

Nm“3 

body  force  on  liquid  water 

Wm-2 

flux  of  heat 

(^h)n'  (^h)n+l 

Wm'^ 

flux  of  heat  within  layers  n  and  n+1  in  the  fringe 

(%)b/  (%)f 

W  m“^ 

flux  of  heat  through  bottom  of  fringe  and  upper 
boundary  zone 

(^h)g/  (^h)w'  (^h)i 

W  K-1  m"l 

thermal  conductivities  of  grains,  water  and  ice 

0 

°C 

temperature 

^0 

K 

absolute  temperature  of  a  flat  ice/water  interface 
in  equilibrium  at  standard  conditions 

h 

Jin-3 

volumetric  latent  heat  of  fusion 

Y 

— 

specific  gravity  of  ice 

z 

m 

position  positive  up  from  the  bottom  of  the  fringe 

{dd/dz%,  (d0/dz)f 

°Cm-l 

temperature  gradient  at  base  of  fringe  and  in  up- 

per  boundary  zone 
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RIGIDICE  Model  of  Secondary  Frost  Heave 

PATRICK  B.  BLACK 


INTRODUCTION 

In  an  earlier  paper.  Black  and  Miller  (1985)  re¬ 
ported  on  progress  to  numerically  solve  the  dif¬ 
ferential  equations  of  secondary  frost  heave  (Miller 
1978)  using  a  simplified  approach.  They  presented 
a  series  of  scenarios  that  predicted  anticipated  re¬ 
sults.  The  program  was  experimentally  confirmed 
by  comparing  observed  behavior  of  independently 
conducted  frost  heave  tests  to  model  predictions 
using  hydrauhc  and  thermal  properties  for  the  test 
soil.  They  found  that  the  model  indicated  that  the 
experimental  data  were  flawed  in  their  stress  mea¬ 
surements,  but  agreed  within  acceptable  tolerance 
with  the  thermal  measurements.  No  additional  ef¬ 
fort  was  made  at  that  time  to  extend  the  analysis 
or  increase  the  utility  of  the  program. 

This  report  is  an  extension  of  the  earlier  work. 
It  presents  improvements  in  two  areas.  First,  a 
minor  mistake  in  the  calculation  of  the  water  flux 
is  corrected;  also,  the  current  formulations  for  the 
scaled  equations  of  ground  freezing  are  used 
(Miller  1990).  Second,  the  computer  code  is  pre¬ 
sented  in  a  form  readily  available  for  enhance¬ 
ments.  This  is  accomplished  by  employing  the 
concepts  of  Object  Oriented  Numerics  (OON), 
which  allow  the  code  to  be  easily  added  to  with¬ 
out  directly  changing  the  original  source  code 
(Wong  et  al.  1993). 

The  concepts  and  equations  of  the  RIGIDICE 
(Black  and  Miller  1985)  model  are  first  reviewed. 
The  new  C-I-+  formulation  of  the  code  is  then  dis¬ 
cussed  and  presented  in  its  entirety  in  Appendix 
A.  The  code  is  linked  as  a  Dynamically  Linked 
Library  (DLL)  that  is  attached  to  MathCad  5.0+ 
(MathSoft  1994).  A  preliminary  sensitivity  study 
is  conducted  to  examine  the  influence  of  uncer¬ 
tainty  in  input  parameters  on  the  variability  of  the 
calculated  output  parameters. 


RIGIDICE 

Black  and  Miller  (1985)  simplified  the  solution 
of  the  differential  equations  for  secondary  frost 
heave  (Miller  1978)  by  employing  two  physical 
assumptions  and  one  numerical  trick.  Their  model, 
as  well  as  this  model,  is  valid  for  one-dimensional, 
incompressible,  air-free  and  solute-free  soil.  The 
solute-free  restriction  assumes  that  no  additional 
osmotic  gradient  is  imposed  over  the  ever  present 
osmotic  gradient  in  the  double  layer  of  the  unfro¬ 
zen  water  surroimding  the  grains.  The  incompress¬ 
ible  restriction  excludes  the  process  of  consolida¬ 
tion  at  this  time.  Finally,  the  air-free  condition  states 
that  the  region  undergoing  heave  must  be  water 
saturated.  This  makes  perfect  sense,  since  air  does 
not  freeze,  which  means  that  the  soil  must  have 
been  saturated,  or  within  10%,  because  of  the  vol¬ 
ume  expansion  necessary  for  a  lens  to  appear. 
These  restrictions  are  minor  since  they  result  in 
the  model  predicting  the  behavior  of  the  highly 
frost-susceptible  fine-grain  silts. 

First  physical  approximation 
(lensing  cycle) 

First,  the  progression  of  frost  heaving  through 
soil  is  approximated  physically  by  a  series  of  in¬ 
dependent  lensing  cycles.  A  lensing  cycle  is  de¬ 
fined  to  be  the  time  step  that  begins  with  the  for¬ 
mation  of  a  lens  and  ceases  with  the  initiation  of  a 
new  lens.  During  each  individual  lensing  cycle, 
the  frozen  soil  above  the  fringe  is  composed  of  a 
series  of  identical  layers  of  ice  and  frozen  soil  as 
depicted  in  Figure  1.  The  thickness  of  these  lenses 
and  the  frozen  soil  between  them  is  equal  to  the 
thickness  of  the  mature  lens  and  interspaced  fro¬ 
zen  soil  at  the  end  of  the  lensing  cycle. 

This  approximation  bounds  the  fringe  between 
the  lower  freezing  front  and  the  upper  boundary 
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zone.  The  freezing  front  is  the  plane  that  separates 
soil  containing  ice  from  the  unfrozen  soil.  The 
upper  boundary  zone  is  not  defined  by  a  plane, 
but  is  rather  a  zone  that  always  contains  a  thick¬ 
ness  of  pure  ice  equal  to  a  mature  ice  lens  and  a 
thickness  of  frozen  soil  equal  to  the  thickness  of 
the  mature  interspace  frozen  soil.  Figure  1  also 
shows  the  location  of  these  regions  at  various 
stages  throughout  the  lensing  cycle. 

Under  steady-state  conditions,  the  inputs  across 
the  freezing  front  are  given  by 


grains 

water 

ice 

sensible  heat 
latent  heat 


GUb 

^sat^b 

0 


and  the  outputs  across  the  boundary  zone  are 


grains 

water 

ice 

sensible  heat 
latent  heat 


Gub 

y^dVb 

/fl^b  + 
(?h)f 

h  [WfV^]  . 


The  global  mass  balance  is  therefore 


('7w)b  +[(^-l)(Wsat -W()]vb  (1) 

and  heat  balance  is 


Mf={cih\  +h[{W,^t-Wf)vi,  +((jw)b]-  (2) 

In  its  purest  sense,  this  results  in  six  unknowns 
and  two  equations.  In  practice,  though,  several  of 
the  variables  are  known  from  other  relationships 
or  observations.  The  numerical  trick  is  to  specify 
at  least  three  of  the  variables  and  guess  another  to 
start  the  numerical  solution  process. 

Temperature  profiles  during  heaving  are  the 
most  common  piece  of  experimental  information 
collected  because  it  is  easy  to  measure.  Freezing 
front  penetration  rate  is  readily  determined  from 
these  data,  so  is  one  logical  variable  to  specify. 
Heat  fluxes  are  also  readily  determined  from  tem¬ 
perature  profile  data  if  the  thermal  conductivities 
are  known.  The  thermal  conductivity  of  the  un¬ 
frozen  soil  is  a  constant  in  this  one-dimensional 
case  because  of  the  air-free  restriction.  The  sen¬ 
sible  heat  flux  across  the  frozen  fringe  there¬ 
fore  makes  another  reasonable  choice.  Unfortu¬ 
nately,  there  is  no  a  priori  justification  for  choos¬ 
ing  from  the  remaining  variables.  The  thermal  con¬ 
ductivity  in  the  frozen  soil  could  be  computed  if 
the  lens  thickness  and  spacing  and  the  residual 
unfrozen  water  content  Wf  were  known.  The  flux 
of  water  into  the  fringe  could  also  be  computed  if 
the  amount  of  heave  were  known.  For  lack  of  any 
further  justification,  the  rate  of  heave  is  fixed. 

If  the  residual  unfrozen  water  Wf  in  the  frozen 
soil  between  ice  lenses  is  known,  then  the  remain- 
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Figure  1.  Lensing  cycle. 
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ing  variables  are  known  for  a  given  soil.  As  a  first 
approximation,  the  amount  of  residual  unfrozen 
water  is  assumed  to  be  equal  to  the  lower  limit  of 
freezing  for  the  soil.  In  some  situations,  this  might 
represent  all  the  required  information  for  the  mod¬ 
eler.  More  often,  the  heaving  pressure,  and  lens 
thickness  and  spacing  are  also  required. 

Second  physical  approximation 
(equivalence  of  instantaneous 
and  averaged  fluxes) 

To  calculate  the  heaving  pressure  and  lens  lo¬ 
cation  within  the  fringe,  profiles  of  temperature 
and  pressures  must  be  calculated  throughout  the 
fringe.  To  perform  such  calculations,  a  method  of 
executing  mass  and  energy  balances  within  the 
fringe  must  be  obtained.  This  is  accomplished  with 
the  second  physical  approximation.  It  states  that 
the  instantaneous  fluxes  of  matter  and  energy  at 
the  beginning  and  the  end  of  the  lensmg  cycle  are 
equal  to  the  averaged  values  of  the  fluxes  during 
the  lensing  cycle. 

Another  way  of  stating  the  second  approxima¬ 
tion  is  that  any  instantaneous  fluctuations  in  mag¬ 
nitude  of  the  mass  and  energy  fluxes  during  the 
lensing  cycle  are  negligible.  This  means  that  the 
magnitude  of  the  penetration  rate,  heave  rate  and 
temperature  gradient  within  the  unfrozen  soil  are 
invariant  with  time  and  space.  In  finite  difference 
form,  the  local  mass  balance  within  the  fringe  is 

(?w)n+i  ”(^w)n 


-Y{vi,+Vi)[W^-W^^r]  (3) 


This  relationship  introduces  the  soil  function  k^, 
the  hydraulic  conductivity.  Each  soil  will  have  its 
unique  hydraulic  conductivity  function.  It  is  not 
a  specific  property  of  heaving  soils,  but  is  rather  a 
general  material  property  of  all  frozen  soils,  heav¬ 
ing  or  non-heaving.  It  must  be  known  in  order  to 
calculate  the  heaving  process. 

To  date,  most  efforts  to  calculate  this  material 
property  are  by  inference.  It  is  inferred  from  ther¬ 
mal  analysis  (van  Loon  et  al.  1988),  back  calcula¬ 
tions  (Ratkje  et  al.  1982  )  and  unsubstantiated  in¬ 
ference  to  non-frozen  soil  (Guymon  and  Luthin 
1974).  Black  and  Miller  (1990)  were  able  to  directly 
measure  the  change  in  hydraulic  conductivity  in 
air-free,  lens-free  and  solute-free  frozen  soil  as  a 
function  of  unfrozen  water  content.  Their  analy¬ 
sis  found  that  if  the  measured  hydraulic  conduc¬ 
tivity  was  expressed  a  function  of  the  difference 
between  the  ice  and  water  pressures 

(t)  =  Wi  -  (6) 

then  the  analogous  expression  given  by  Brooks 
and  Corey  (1964)  for  partially  saturated  and  ice- 
free  soil  could  be  transformed  to 


and 


W(^)  =  (W3at-Wd) 


u 


Y 

) 


+  Wd 


(^w)sat 


^<l>b 


(7) 

(8) 


Fourier's  law 

The  flux  of  thermal  energy  through  the  fringe 
is  assumed  to  follow  the  Fourier's  law 


and  thermal  balance  is 

The  remaining  information  required  to  com¬ 
plete  all  calculations  are  statements  for  water  and 
ice  pressures,  temperature  and  a  criterion  for  lens 
initiation. 


—  =  -51l  (9) 

dz  fch‘ 

This  relationship  introduces  another  soil  function, 
fch,  the  thermal  conductivity  for  the  frozen  soil. 
Just  as  the  hydraulic  conductivity,  this  material 
property  is  also  a  function  of  the  pressure  differ¬ 
ence  between  the  ice  and  water  pressures.  One 
standard  expression  is  the  geometric  mean  formu¬ 
lation  of  Farouki  (1981) 


Darcy's  law 

The  flux  of  water  through  the  fringe  is  assumed 
to  obey  Darcy's  law 


dWw 

dz 


“/w 


(5) 


Clapeyron  equation 

The  equilibrium  condition  between  the  ice  and 
water  pressures  and  temperature  is  given  by  the 
Clapeyron  equation 
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dz  ^  ’  dz  -Go  dz'  ^  ’ 

Criterion  for 
lens  initiation 

To  calculate  when  a  new  lens  forms  and  the 
lensing  cycle  is  complete,  a  criterion  for  lens  ini¬ 
tiation  is  required.  Simply  stated,  a  new  lens  will 
form  when  ice  can  penetrate  between  the  grains, 
causing  the  grains  to  separate.  Other  factors  must 
be  included  if  the  porous  material  was  not  granu¬ 
lar.  For  example,  cements  and  rocks  that  have 
chemical  bonds  between  the  grains  must  have  the 
bonds  broken  before  the  ice  can  penetrate. 

Geotechnical  engineers  express  the  state  of 
stress  within  a  granular  material  by  the  Terzaghi 
equation 


aT=Oe+an  (12) 

and  the  condition  when  the  grains  separate  by 
CTe=0.  (13) 

At  this  instance 


01=  On-  (14) 

Now  if  CTj  is  assumed  to  be  equal  to  the  maxi¬ 
mum  heaving  pressure,  then  the  location  where  a 
new  ice  lens  forms  is  where  also  equals  the 
maximum  heaving  pressure.  Snyder  and  Miller 
(1985)  found  that  they  could  correctly  model  their 
empirically  measured  neutral  stress  data  by 

On=XWw+(l-xK  (15) 

in  which  the  new  soil  function  was  determined  to 
be 


X  ((!>)== 


1 

2 


W(4>)-Wd 

IVsat-Wd 


0.3  «  W(<l>j)-Wd 
h  ^sat-Wd 


(16) 


Algorithm 

The  solution  strategy  begins  by  stating  values 
for  Pj,  and  (d^ldz)^  and  all  necessary  param- 
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Figure  2.  Flow  chart  ofRIGIDICE. 
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eters  for  the  soil  functions.  Next,  the  residual  un¬ 
frozen  water  content  in  the  mature  frozen  soil  is 
assumed  to  be  equal  to  the  lower  limit  of  freez¬ 
ing.  Equations  1  and  2  are  then  solved  for  the  water 
flux  into  the  fringe  and  the  flux  of  heat  out  of  the 
upper  boundary. 

Profiles  of  pressures,  stresses  and  temperature 
are  then  calculated  within  the  fringe  starting  at 
the  freezing  front.  The  value  of  (j)^  at  the  freezing 
front  is  obtained  by  setting  the  local  water  pres¬ 
sure  to  zero.  This  is  equivalent  to  assuming  that 
the  water  table  is  at  the  freezing  front.  If  the  water 
pressure  is  not  zero,  then  the  resulting  heaving 
pressure  obtained  by  assuming  zero  must  be  ad¬ 
justed  using  eq  11. 

Calculations  are  done  in  terms  of  (|).  In  finite 
difference  form,  balances  are  conducted  across  lay¬ 
ers  of  constant  d(j).  This  has  the  benefit  of  generat¬ 
ing  thin  spatial  layers  where  (j)  is  changing  rap¬ 
idly  and  large  spatial  layers  where  it  changes  least. 

Darcy's  law  gives  the  value  of  the  water  pres¬ 
sure  across  the  layer  and  Fourier's  law  gives  the 
temperature.  The  Clapeyron  equation  then  gives 
the  ice  pressure  as  well  as  the  thickness  of  the  layer. 
Neutral  stress  is  calculated  for  each  layer  and  a 
running  account  of  its  magnitude  recorded  to  de¬ 
termine  its  maximum  value.  A  running  sum  of  the 
unfrozen  water  content  is  also  made,  starting  from 
the  location  of  maximum  neutral  stress  to  deter¬ 
mine  the  residual  water  content.  Calculations  con¬ 
tinue  until  the  ice  pressure  equals  the  maximum 
neutral  stress.  This  is  the  location  of  the  base  of 
the  lens  and  the  new  lens  will  form  where  the  neu¬ 
tral  stress  was  maximum.  The  distance  between 
these  two  locations  gives  the  lens  spacing.  The  cal¬ 
culated  residual  water  content  is  then  compared 
to  the  original  guess.  If  the  difference  is  unaccept¬ 
able,  the  layer  by  layer  calculations  are  performed 
again  with  the  new  guess  until  the  resulting  change 
in  residual  water  content  is  acceptable. 

A  flow  chart  representation  of  this  algorithm 
is  presented  in  Figure  2. 

RIGIDICE  WITH 

OBJECT-ORIENTED  NUMERICS  (OON) 

Code  was  written  to  solve  eq  1  through  16  us¬ 
ing  the  algorithm  outlined  above.  To  allow  for  the 
greatest  flexibility  of  use,  as  well  as  ease  of  future 
enhancements,  C++  was  used.  The  entire  code  is 
presented  in  Appendix  A. 

The  benefit  of  using  OON  is  that  it  allows  the 
writing  of  code  in  separate  layers  that  are  easily 


merged  together.  This  might  be  interpreted  as  the 
standard  approach  of  writing  a  series  of  subrou¬ 
tines,  but  it  is  more.  The  traditional  approach  is  to 
write  a  series  of  procedures  to  numerically  solve 
the  problem.  These  procedures  are  sequentially 
solved  and  in  large  programs  result  in  many  pages 
of  critically  linked  lines  of  code.  When,  for  ex¬ 
ample,  the  routine  to  calculate  hydraulic  conduc¬ 
tivity  must  be  changed,  that  section  of  code  must 
be  found  and  modified.  Care  must  be  taken  not  to 
remove  variables  that  are  used  by  the  rest  of  the 
program  as  well  as  not  to  add  variables  that  are 
already  being  used  in  other  parts  of  the  program. 

The  OON  approach  is  to  break  the  problem 
down  into  separate  self-contained  units  called 
classes,  a  class  being  a  collection  of  data  and  op¬ 
erations.  Data  in  one  class  can  be  made  inacces¬ 
sible  to  other  classes  to  prevent  inadvertent 
changes  by  later  modifications  to  the  program. 
Additional  classes  can  be  made  that  inherit  the 
properties  of  existing  classes.  Operations  used  by 
a  class  can  be  expanded  by  a  process  called  poly¬ 
morphism.  Again,  if  a  new  function  to  calculate 
hydraulic  conductivity  is  needed,  it  is  not  neces¬ 
sary  to  write  a  completely  new  class,  but  just 
merely  to  add  to  the  current.  In  other  words,  the 
new  function  would  access  all  the  data  used  by 
the  old  as  well  as  any  new  data  that  it  requires, 
and  the  new  data  are  prevented  from  interfering 
with  any  existing  data  in  the  program.  The  tradi¬ 
tional  approach  leads  to  mistakes  and  debugging 
problems.  The  OON  approach  rests  on  the  belief 
that  if  the  original  code  worked,  then  don't  change 
it,  just  add  improvements. 

RIGIDICE  is  setup  as  a  series  of  C++  classes. 
At  the  lowest  level,  functions  and  data  associated 
with  a  particular  soil  are  grouped  together  in  a 
class  called  Tsoil.  Functions  and  data  associated 
with  the  boundary  conditions  are  grouped  together 
as  a  class  Tbnds  and  numerical  precision  in  a  class 
Ttol  Since  all  variables  are  reduced  to  dimension¬ 
less  form  to  help  in  numerical  calculations,  a  class 
that  does  scaling  is  called  REDUCE.  Finally,  the 
calculations  that  perform  the  necessary  algorithm 
are  made  in  the  class  Trigidice.  This  class  is  de¬ 
clared  so  that  it  inherits  all  the  properties  of  the 
other  classes.  If  a  different  algorithm  is  necessary, 
then  a  new  one  can  be  written  that  also  inherits 
the  other  classes  while  never  touching  the  origi¬ 
nal  source  code  for  the  other  primitive  classes. 

In  total,  the  program  requires  the  initiation  of 
18  variables.  It  then  returns  the  values  of  eight 
calculated  variables.  Figure  3  shows  the  input 
screen  of  the  MathCad  program  that  runs  RIGIDICE. 
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This  program  calculates  the  rate  of  heave  for  a  given  set  of  initial  and  boundary  conditions.  The  strategy  is  the  same  as 
presented  by  Black  and  Miller  (1 985).  It  Is  assumed  that  the  heaving  process  can  be  approximated  with  a  time  step  called 
the  lensing  cycle  and  that  any  fluctuations  during  this  lensing  cycle  average  out  to  values  that  represent  the  true  behavior 
during  the  lensing  cycle. 

Scaling  factors 

fw:  =  1 

Thermal  conductivities 

Khw:  =  0.52 
Khi:  =  2.32 
Khg:  =  3.42 

Soil  water  properties 

Wsat:  =  0.42 
Wd:  -  0.02 
Kwsat:  =  1.0  •  10"® 

Brooks  and  Corey  functions  settings  (heave,  penetration):  = 

(J)b:=  11.196 
a:  =  0.36 
P:  =  2.6 

Boundary  conditions 

heave:  =  10 
penetration:  =  100 
gradOunfrozen:  =  -10 

Numeric  settings 

prec:  =  0.1 
resol:  =  0.01 
MaxLayers:  =  1 00000 

The  function  RIGIDICE  takes  the  contents  of  settings  and  returns  a  1-dimensional  array  containing  (heave  pressure, 
temperature  gradient  in  frozen  soil,  heat  fluxes  into  the  fringe  and  out  of  the  upper  boundary  zone,  the  water  flux  into 
the  fringe  and  the  lens  spacing  and  thickness  and  iterations) 

I 

I  (ui  gradOf  qhin  qhout  qwin  spacing  thickness  iters):  =  RIGIDICE  (settings  [heave,  penetration]) 

jui  =  76.06  spacing  =  0.118  thickness  =  4.251  qhin=  15.504  qhout  =  199.533  iters  =  2  gradOf  = -70.255 

Figure  3.  MathCad  program  to  run  RIGIDICE, 


X 

fw 

Khw 

Khi 

Khg 

Wsat 

Wd 

Kwsat 

a 

P 

({)b 

heave 

penetration 

gradOunfrozen 

prec 

resol 

_  MaxLayers 
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The  user  is  free  to  modify  the  contents  of  all  18 
input  parameters  that  are  passed  to  the  array  set¬ 
tings.  The  function  RIGIDICE  (settings)  returns  an 
array  that  contains  the  eight  calculated  variables. 


DISCUSSION 

Owing  to  the  number  of  input  parameters,  a 
thorough  sensitivity  study  of  the  program  is  very 
time  consuming.  This  study  employed  a  quicker 
approach  in  which  the  relative  significance  of  the 
various  input  parameters  is  obtained  through  a 
simple  scheme  of  sequentially  varying  one  param¬ 
eter  at  a  time  and  noting  the  resulting  effect  on 
the  important  calculated  parameters.  Each  input 
parameter  is  modified  by  10%  and  the  calculated 
results  noted.  The  initial  reference  value  for  each 
parameter  is  obtained  from  past  measurements  for 
one  particular  silty  soil  (Black  and  Miller  1985, 1990). 

It  is  evident  that  the  heaving  pressure  is  of  pri¬ 
mary  concern  to  the  end  user.  The  heaving  pres¬ 
sure  is  therefore  the  output  parameter  to  be  ex¬ 
amined.  In  addition  to  examining  how  the  heave 
pressure  alone  changes  with  each  parameter,  the 
heave  rate  and  heave  pressure  behavior  will  be 
examined  as  the  other  parameters  are  modified 
one  at  a  time.  This  approach  is  similar  to  how  an 


end  user  might  use  the  program.  By  stating  all 
the  other  parameters,  the  end  user  is  then  able  to 
predict  heave  rate  as  a  function  of  heave  pressure 
(i.e.,  overburden  pressure). 


0  20  40  60  80  100 

Vi,  Heave  (mm/day) 

c.  Saturated  hydraulic  conductivity,  (Rjo^sat- 


Figure  4.  Changes  in  heave  pressure  vs.  heave  rate  for  changes  in  chosen  input  parameters. 
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Figure  4  (coni' d).  Changes  in  heave  pressure  vs. 

The  strategy  is  to  examine  the  influence  that 
changes  in  the  12  chosen  input  parameters  have  on 
the  magnitude  of  the  predicted  heaving  pressure. 
This  is  accomplished  by  running  simulations  that 
vary  the  value  of  these  parameters  by  10%  from  the 
initial  reference  values  in  Table  1.  The  resulting  pre- 


g.  Penetration  rate,  v^. 

ve  rate  for  changes  in  chosen  input  parameters. 

dieted  heaving  pressures  are  then  computed  and 
displayed  in  Table  1  for  direct  comparison  along  with 
a  series  of  graphs  that  show  the  overall  effect  on  the 
heave  pressure  as  a  function  of  heave  rate. 

This  approach  is  important  in  designing  experi¬ 
mental  tests  to  aid  in  the  development  of  engi- 
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Figure  4  (cont'd). 


neering  design  criteria.  Variables  that  cause  the 
greatest  change  in  predicted  heaving  pressure  are 
assumed  to  be  important.  The  end  user  should 
therefore  concentrate  on  obtaining  correct  values 
for  these  important  variables  if  the  predicted  re¬ 
sults  are  to  be  useful. 


Figures  4a,  b  and  c  show  how  a  10%  change  in 
either  direction  from  the  reference  value  for  the 
three  soil  water  parameters  influences  the  calcu¬ 
lated  heave  pressure  for  a  given  heave  rate.  In  ei¬ 
ther  case,  the  resulting  deviation  is  less  than  the 
10%  change  in  the  control  parameter.  A 10%  change 
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Table  1.  Sensitivity  results  of  calculated  heave  pressure  to  input  parameter  changes  in  RIGIDICE. 


_ Calculated  heave  pressure  (kPa) 

0.9  Reference  1.2  Reference 

Program  parameter _ Reference  value _ relative  error*  Reference  relative  error* 

Soil  water  properties 


Wsat 

0.42  (m^  m“^) 

70.22 

-7.7 

76.06 

82.05 

7.9 

0.02  (m^  m~^) 

76.08 

2.6 

76.06 

75.72 

0.4 

(^w)sat 

1x10-8  (m  s-i) 

72.33 

-4.9 

76.06 

79.58 

4.6 

Brooks  and  Corey  constants 

11.196  (kPa) 

67.94 

-10.7 

76.06 

81.82 

7.6 

a 

0.36 

71.55 

-5.9 

76.06 

79.86 

5.0 

p 

2.6 

109.03 

43.3 

76.06 

56.86 

-25.2 

Boundary  conditions 

10  (mm  day“^) 

81.38 

7.0 

76.06 

71.39 

-6.1 

100  (mm  day“^) 

71.59 

-5.9 

76.06 

80.24 

5.5 

V0 

-10  (°C  m-i) 

75.47 

-0.8 

76.06 

76.38 

0.4 

Thermal  conductivities 

0.52  (W/K-l  m-l) 

76.38 

0.4 

76.06 

75.38 

-0.9 

<Kh 

2.32  (W/K-l  m-l) 

76.74 

0.9 

76.06 

75.12 

-1.2 

3.42  (W/K-l  m-i) 

77.86 

2.4 

76.06 

74.33 

-2.3 

Scaling  factors 

X 

1.0x10-^  (m) 

c 

1.0x10-2  (m) 

/w 

1  (m  s-2) 

Numeric  settings 

Precision 

0.1 

Resolution 

0.01 

Maximum  layers 

1.0x10^ 

*  ,  .  /  N  value  -  reference 

relative  error  [%)  - - - - 100 

reference 
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in  the  saturated  water  content  Wgat  resulted  in  the 
largest  variation  of  the  three,  as  listed  in  Table  1, 
with  a  relative  error  range  of  -8  to  8%.  These  re¬ 
sults  suggest  that  a  10%  uncertainty  in  the  magni¬ 
tude  of  Wsat/  VVd  ai'd  (*^w)sat  should  not  have  a 
severe  adverse  effect  on  the  predicted  heave  pres¬ 
sure. 

The  largest  relative  errors  result  from  a  10% 
change  in  the  Brooks  and  Corey  constants  (j)b  and 
P  for  the  soil  as  shown  in  Figures  4d  and  4f.  Table 
1  lists  the  largest  range  of  43  to  -25%  relative  er¬ 
ror  in  the  calculated  heave  pressure  for  a  10% 
change  in  P  and  a  —10  to  8%  error  for  (j)b-  Clearly, 
this  indicates  that  the  hydrologic  properties  of  the 
soil  must  be  known  if  any  accurate  heave  behav¬ 
ior  is  to  be  calculated.  This  is  unfortunate  since 
this  information  is  rarely  even  collected.  Figure 
4e  and  Table  1  shows  a  smaller  response,  similar 
to  the  soil  water  properties,  for  a  change  in  a. 

Surprisingly,  a  10%  change  in  the  boundary 
conditions  also  gives  a  similar  response  in  rela¬ 
tive  error  to  heave  pressure  as  the  soil  water  prop¬ 
erties.  Figures  4g  and  h  show  the  response  for  a 
change  in  penetration  V\)  and  the  temperature  gra¬ 
dient  in  the  unfrozen  soil  V0.  The  heave  rate 
response  to  calculated  heave  pressure  Wj  is  obtained 
from  any  of  the  graphs  as  it  is  always  the  center 
reference  line.  Again,  this  is  unfortunate  since  the 
penetration  rate  and  temperature  gradient  are  eas¬ 
ily  obtained  from  temperature  profiles.  This  means 
that  more  accurate  temperature  measurements  will 
not  necessarily  result  in  better  predictive  capabil¬ 
ity. 

The  last  group  of  graphs  displays  the  response 
to  an  imcertainty  in  thermal  conductivities.  Inspec¬ 
tion  of  Figures  4i,  j  and  k  along  with  Table  1  re¬ 
veals  that  an  imcertainty  in  this  group  of  param¬ 
eters  has  the  least  influence  of  all  parameters  on 
the  calculated  heave  pressure. 


CONCLUSIONS 

The  past  problem  of  all  frost  heave  models  was 
the  trade-off  between  physical  correctness  and  ease 
of  use.  Those  models  that  are  easy  to  calculate  tend 
to  be  based  upon  curve  fitting  heave  experiments 
(Blanchard  and  Fremond  1985)  or  incorrectly  ap¬ 
plying  phase  equilibrium  (Harlan  1973).  The  mod¬ 
els  that  are  physically  based  are  difficult  to  imple¬ 
ment  and  require  time-consuming  computation 
(O'Neill  and  Miller  1985).  Past  efforts  to  overcome 
this  dilemma  relied  on  modifying  the  original 
equations  of  Miller  by  simplifying  the  physics  to 


make  the  mathematics  trivial  (Gilpin  1980,  Holden 
1983)  or  attempting  to  maintain  the  physics  while 
leading  to  a  reasonable  solution  strategy  (Black 
and  Miller  1985).  This  later  approach  was  success¬ 
fully  employed  in  this  paper. 

The  differential  equations  of  secondary  frost 
heave  (Miller  1978)  are  numerically  solved  in  fi¬ 
nite  difference  form  with  the  program  RIGIDICE. 
By  choosing  the  language  C++  and  its  object  ori¬ 
ented  nature,  the  core  program  is  easily  attached 
to  the  mathematical  analysis  program  MathCad 
5.0+.  The  ease  of  analysis  in  this  new  format  al¬ 
lows  simulations  to  be  conducted  in  any  physi¬ 
cally  correct  model  that  the  end  user  requires. 

A  brief  sensitivity  analysis  of  several  impor¬ 
tant  parameters  shows  that  current  practices  of 
ground  freezing  monitoring  are  flawed.  Calculated 
behavior  for  heave  rate  and  pressure  was  found 
to  be  largely  insensitive  to  penetration  rate  and 
temperature  gradients.  Likewise,  unfrozen  water 
content  uncertainty  did  not  strongly  influence  the 
heave  rate  and  pressure  behavior.  Hydraulic  con¬ 
ductivity,  though,  was  found  to  have  a  dramatic 
influence. 

These  results  indicate  that  the  effort  expended 
in  making  accurate  temperature  profiles  in  freez¬ 
ing  ground  is  not  necessary.  It  also  indicates  that 
the  great  efforts  required  to  monitor  unfrozen  wa¬ 
ter  content  changes  in  the  field  might  also  be  un¬ 
necessary.  The  sensitivity  study  did  demonstrate 
the  importance  of  the  hydraulic  conductivity.  A 
wise  use  of  research  and  development  efforts 
should  therefore  be  to  develop  techniques  to  mea¬ 
sure  the  hydraulic  conductivity  in  the  laboratory 
and  monitor  it  in  the  field. 

The  current  model  relies  upon  assumed  stan¬ 
dard  expressions  for  thermal,  hydraulic  and  stress 
behavior  in  soil.  Other  expressions  need  to  be  ex¬ 
plored  (i.e.,  van  Genuchten's  expression  for  hy¬ 
draulic  properties)  in  order  to  generalize  the  ap¬ 
plicability  of  this  model.  The  OON  approach  em¬ 
ployed  will  make  this  a  reasonable  task.  Empiri¬ 
cal  testing  of  this  model  is  currently  under  way 
with  a  refrigerated  centrifuge.  This  approach  will 
test  the  scaling  laws  for  freezing  (Miller  1990)  as 
well  as  the  predictions  of  this  model. 
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APPENDIX  A:  RIGIDICE  MODEL 


RI6IDICE.H 


0 

1 

2 

3 
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//  acceleration  due  to  gravity 

#define  grav  ((const  double)  9.800)  //  ni/sec**2 

//  density  of  liquid  water 

#define  rhoW  ((const  double)  1.000E03)  //  kg/m**3 

//  density  of  ice 

#define  rhol  ((const  double)  0.917E03)  //  kg/rti**3 

//  viscosity  of  water 

#define  nu  ((const  double)  1.787E-3)  //  kg/ni*s 

//  specific  surface  free  energy  of  a  ice-air  interface 
0  #define  gammalA  ((const  double)  0.100)  //  n/m 

//  Latent  heat  of  fusion  per  unit  volume  of  melt 
#define  hlW  ((const  double)  3.335E05)  //  j/m**3 

//  standard  melting  point  of  bulk  ice,  or  its  analogue,  when  exposed 

//to  air  at  standard  atmospheric  pressure 

#define  thetaO  ((const  double)  273.15)  //  kelvins 

//  conversions 

//  mm/day  /  (m/s) 

#define  mmpdTomps  ((const  double)  8.64e7) 

//  kPa  to  Pa 

#define  kP2P  ((const  doiible)  l.OOeS) 


class  REDUCE! 
protected: 

double  orglamda, 

orgeta, 
orgFW; 

public: 

//  Constructor 

REDUCE ( ) ; 

void  InitScales ( ) ; 

void  FirstScales (double  micro,  double  macro,  double  body) 
double  lamda; 

double  eta; 

//  double  body; 

double  FW; 
double  Y; 

double  L (double  rSI) ; 
double  P (double  pSI)  ; 

//  reduced  latent  heat  of  fusion,  per  unit  volume  of  melt 
double  H; 

double  RTHETA( double  thetaSI) ; 
double  SIGMA(double  sigmaSI) ; 
double  GRAD (double  gradSI) ; 
double  F(double  fSI) ; 
double  KW( double  kwSI) ; 
double  V( double  vSI) ; 
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double  DIV{double  divSI) 
double  T{double  tSI) ; 
double  C{double  cSI)  ; 
double  KH(double  khSI) ; 
double  QH( double  qhSI) ; 


class  Tsoil{ 
public : 


REDUCE  R; 

//  saturated  volumetric  water  content 
double  WSAT; 

//  lower  limit  of  drying 
double  WD; 

//  reduced  saturated  hydraulic  conductivity 
double  KWSAT; 

public : 

//  Constructor 

Tsoil ( ) ; 

void  In itWATER (double  wsat,  double  wd,  double  ksat) 
//  calculates  water  content  from  degree  of  saturation 
double  W  (double  DSAT) ; 

private : 

//  volume  fraction  of  ice 
double  I; 

//  volume  fraction  of  soil 
double  G; 

//  thermal  conductivity  of  water 
double  KHW; 

//  thermal  conductivity  of  ice 
double  KHI; 

//  thermal  conductivity  of  soil  grains 
double  KHG; 

public : 

void  Ini tTherm( double  khw,  double  khi,  double  khg) ; 
//  calculates  thermal  conductivity  for  given  composition 
double  RKH  (double  W) ; 

//  returns  the  protected  reduced  ice  conductivity 
double  RKHI ( ) ; 


protected: 
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//  unfrozen  water  content  exponent 
double  ALPHA; 

private: 

//  frozen  hydraulic  conductivity  exponent 
double  BETA; 

protected: 

//  reduced  ice  intrusion  pressure 
double  PHIB; 

public : 

void  InitBC (double  alpha,  double  beta,  double  phib) ; 
double  DSAT (double  PHI) ; 
double  RKW(double  PHI) ; 

/ /  running  sum 
private : 

double  PROD; 
int  FlagSny; 

public : 

//  initialize  starting  parameters 
void  InitStress ( )  ; 

//  calculates  Snyder's  factor 

double  SNYDER  (double  PHI,  double  DSAT,  double  deltaDSAT) ; 

}; 


//  boundary  conditions 
class  Tbnds 

{ 

public:  REDUCE  R; 
protected: 

/ /  heave  rate 
double  VI,  VIorg, 

//  penetration  rate 
VB,  VBorg, 

//  unfrozen  temperature  gradient 
GRADTB ,  GRADTBorg , 

//  sought  after  heave  pressure 
Ulorg; 


public : 

void  InitQuasi (double  press,  double  pene,  double  gradtb) ; 
void  Ini tContinuum (double  heave,  double  pene,  double  gradtb) 
Tbnds ( ) ; 

); 


//  variables  used  for  numeric  precision 
class  Ttol 

{ 

protected: 

double  orgresolution, 

orgprecision; 

public : 
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double  precision, 
resolution, 

LCMAX; 


//  check  on  Wd 

If  layer  thickness 

//  cutoff  number  of  layers 
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double  DPHI; 
void  InitTol { ) ; 

void  FirstTol (double  prec,  double  resol,  double  Icmax) ; 
Ttol  (  )  ; 

}; 


//  variables  used  in  calculations 
class  Trigidice: 

public  Tbnds, 
public  Ttol, 
public  REDUCE, 
public  Tsoil 
{ 

public : 

double  QWl,  QW2,//  reduced  mass  flux 

DW,//  reduced  change  in  water  content  across  layer 
GRADUWl ,  GRADUW2,//  reduced  water  pressure  gradient 
QHl,  QH2,//  reduced  thermal  flux 

GRADTHETAl,  GRADTHETA2 , / /  reduced  thermal  gradient 
PHI,//  reduced  phi 

GRADPHIl,  GRADPHI2,//  reduced  gradient  of  PHI  from  Clapeyron 

Z,//  reduced  thickness 

ZN, //  reduce  location  of  new  lens 

ZI,//  reduce  location  of  old  lens 

DZ,//  reduced  layer  thickness 

UW, //  reduced  water  pressure 

DUW, //  reduced  change  in  water  pressure  across  layer 

UI,//  reduced  ice  pressure 

CHI ,/ /reduced  stress  partition  factor 

UN,//  reduced  neutral  stress 

UNold,//  reduced  previous  neutral  stress 

UNmax,//  reduced  maximum  neutral  stress 

THETA,//  reduced  temperature 

DTHETA,//  reduced  temperature  change  across  layer 

WF,//  unfrozen  water  in  mature  frozen  zone 

FSI; 

Trigidice { ) ; 

void  LAYER (double  rphi,  double  rdphi) ; 
void  PROFILEO; 
void  GLOBAL ( ) ; 

double  oldPHI; 
double  sum; 
int  iteration, 

//  error  signal  for  GRADPHI  <  0 
GRADPROB, 

//  not  enough  layers 
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LAYERPROB; 
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217 

218 

219 

220 
221 
222 

223  }; 


double  tic, 

//  number  of  layers 
counter; 
double  maybe; 
int  check ( ) ; 
int  prec ( ) ; 
void  errors { ) ; 
void  once ( ) ; 
int  pressdif f ( ) ; 
int  slopeO; 
void  quasi ( ) ; 
void  InitTCALCO; 
double  heave  0 ; 
double  iter out ( ) ; 
double  pressure ( ) ; 
double  tgradout ( ) ; 
double  heat in ( ) ; 
double  heatout ( ) ; 
double  waterinO; 
double  spacing ( ) ; 
double  thickness { ) ; 
void  lookO; 
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#include  <iostream. h> 

#include  <math.h> 

#include  "rigidice.h" 

/  /  REDUCED  VARIABLES 

//  From: 

//  R.  D.  Miller  (1990)  Scaling  of  freezing  phenomena, 

//  in  Scaling  in  soil  physics :principles  and  applications 
//  SSSA  special  Publication  Number  25 
0 

REDUCE : : REDUCE ( )  { 

//  specific  gravity  of  ice? 

Y=rhoI/rhoW; 
orglamda=l . Oe-6 ; 
orgeta=0.01; 

FW  =  1; 

orgFW  =  F{FW* (-rhoW*grav) ) ; 

//  reduced  latent  heat  of  fusion,  per  unit  volume  of  melt 
H  =rhoW* ( lamda/gammalA) *hIW; 

InitScales ( ) ; 

} 

void  REDUCE :: FirstScales {double  micro,  double  macro,  double  body) { 
orglamda  =  micro; 
orgeta  =  macro; 

FW  =  body; 

orgFW  =  F(FW* (-rhoW*grav) ) ; 

//  reduced  latent  heat  of  fusion,  per  unit  volume  of  melt 
H  =rhoW* ( lamda/gammalA) *hIW; 

InitScales ( ) ; 

} 


void  REDUCE :: InitScales () { 
lamda  =  orglamda; 
eta  =  orgeta; 

FW  =  orgFW; 

} 


//  reduced  length 

double  REDUCE: :L (double  rSI) 

{  return  ( ( 1/lamda) *rSI ) ; } 

//  reduced  pressure 

double  REDUCE: :P (double  pSI) 

{  return  ( ( lamda/gammalA) *pSI) ; } 
/ /  reduced  temperature 

double  REDUCE: :RTHETA( double  thetaSI) 

{  return  {( 1/thetaO )* thetaSI ); } 


//  reduced  stress 

double  REDUCE :: SIGMA (double  sigmaSI) 

{  return  ( ( lamda/gammalA) *sigmaSI) ; } 

//  reduced  gradient 
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double  REDUCE: : GRAD (double  gradSI) 

{  return  {eta*gradSI) ; } 

//  reduced  body  force 

double  REDUCE: :F (double  fSI) 

{  return  ( ( lamda* eta /gamma lA) *fSI) ; } 
//  reduced  capillary  conductivity 

double  REDUCE : :KW( double  kwSI) 


{  return  ( (nu/ ( lamda* lamda) ) *kwSI) ; } 

reduced  velocity 

double  REDUCE: :V(double  vSI) 

{  return  ( (eta*nu/ ( lamda*gammaIA) ) *vSI) ; } 

reduced  divergence 

double  REDUCE: :DIV( double  divSI) 

{  return  (eta*divSI ) ; } 

reduced  time 

double  REDUCE: :T (double  tSI) 

{  return  ( (lamda*gammaIA/ (eta*eta*nu) ) *tSI) ; } 
reduced  volumetric  heat  capacity 

double  REDUCE: :C (double  cSI) 

{  return  ( ( lamda* thetaO /gammalA) *cSI) ; } 
reduce  thermal  conductivity 

double  REDUCE: :KH( double  khSI) 

{  return  ( (nu*theta0/ (gammaIA*gammaIA) ) *khSI) ; ) 

reduced  heat  flux 

double  REDUCE: :QH( double  qhSI) 

{  return  (eta*nu/ (gammaIA*garamaIA) *qhSI) ; } 


SOIL  FUNCTIONS 


Tsoil : : Tsoil ( ) { 

InitWATER (0.42,0.02,1. 0e~8 )  ; 
InitBC (0.36,2.6,11.196); 
InitTherm(0 .52,2.32,3.42) ; 
InitStress ( ) ; 

} 


void  Tsoil :: InitWATER (double  wsat,  double  wd,  double  ksat) { 
//  saturated  volumetric  water  content 
//  WSAT  =  0.42;^ 

WSAT  =  wsat; 
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//  lower  limit  of  drying 
//  WD  =  0.02; 

WD  =  wd; 

//  saturated  hydraulic  conductivity 
//  KWSAT=1.0e-8;  //  m/s 

KWSAT  =  ksat; 

KWSAT=KWSAT/ (rhoW*grav) ; 

KWSAT  =  R.  KW  ( KWSAT  ); 

) 

//  given  degree  of  saturation 
double  Tsoil : :W (double  DSAT) { 

return (  DSAT* (WSAT-WD) +WD) ; 

} 

// 

// 

//  Thermal  conductivity 

void  Tsoil :: InitTherm (double  khw,  double  khi,  double  khg) { 

//  volume  fraction  of  soil 
G=1-WSAT; 

//  volume  fraction  of  ice 
I=:1-G-WD; 

//  thermal  conductivity  of  water 
//  KHW  =0.52;  //  W/m 

KHW  =  khw; 

KHW  =  R.KH(KHW) ; 

//  thermal  conductivity  of  ice 
//  KHI  =2.32;  //  W/m 

KHI  =  khi; 

KHI  =  R.KH(KHI); 

//  thermal  conductivity  of  soil  grains 

//  KHG=3.42;  //  W/m 

KHG  =  khg; 

KHG  =  R.KH(KHG); 

} 

//  calculates  thermal  conductivity  for  given  composition 
double  Tsoil: :RKH  (double  W) { 

I  =  1  -  G  -  W; 

return  (  pow(KHG,G) *pow(KHW, W) *pow(KHI, I) ) ; 

} 

//to  get  KHI 
double  Tsoil: :RKHI() { 

return  (KHI) ; 

} 

// 

// 

//  Brooks  &  Corey  equations 

// 

void  Tsoil :: InitBC (double  alpha,  double  beta,  double  phib)  { 
II  ALPHA  =  0.36; 
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//  BETA  =  2.6; 

ALPHA  =  alpha; 

BETA  =  beta; 

//  ice  entry  pressure 
//  PHIB  =  11.196*kP2P;  //  kPa 

PHIB  =  phib; 

PHIB  =  R.P(PHIB*kP2P)  ; 

} 

//  returns  degree  of  saturation  for  a  given  phi 
double  Tsoil: :DSAT{ double  PHI)  { 

return  (pow( (PHIB /PHI) , ALPHA) ) ; 

} 

//  hydraulic  conductivity  for  a  give  phi 
double  Tsoil: :RKW{ double  PHI)  { 

//  return  (be .KWSAT*pow( (PHIB /PHI ), BETA) ) ; 

return  (KWSAT*pow( (PHIB/PHI) .BETA) ) ; } 


Stress  partition  factor 


void  Tsoil :: InitStress 0 { 
set  defaults  to  obtain  a  Snyder  factor  of  1 
PROD=0.0; 

FlagSny  =  0; 


calculates  Snyder's  factor 

double  Tsoil: : SNYDER  (double  PHI, 


double  DSAT, 


double  deltaDSAT) { 

if  (FlagSny  ==  0)  { 

FlagSny  =  1 ; 

PROD  =  0; 
return  ( 1) ; 

} 

else{ 

PROD  =  PROD  +  PHI*deltaDSAT; 
return(  0.5*  (DSAT- (0 . 3/PHI) *PROD) ) ; 

) 

} 


Tbnds : : Tbnds ( ) { 

InitQuasi ( 100 ,10,-10) ; 
InitContinuuin(10 , 100  , -10)  ; 

} 


void  Tbnds :: Ini tContinuum( double  heave,  double  pene,  double  gradtb) { 
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204 

VI  =  heave; 

205 

VIorg  =  R. V (VI/mmpdTomps ) ; 

206 

VI  =  VIorg; 

207 

VB  =  pene; 

208 

VBorg  =  R. V(VB/iTimpdTomps)  ; 

209 

VB  =  VBorg; 

210 

GRADTB  =  gradtb; 

211 

GRADTBorg  =  R . GRAD (R . RTHETA ( gradtb) ) ; 

212 

GRADTB  =  GRADTBorg; 

213 

} 

214 

215 

void  Tbnds 

:: Ini tQuasi (double  press,  double  pene,  double  gradtb)' 

216 

217 

Ulorg  =  R. P (press*kP2P) ; 

218 

VBorg  =  R.  V  (pene/nunpdTomps  )  ; 

219 

GRADTBorg  =  R . GRAD (R . RTHETA ( gradtb) ) ; 

220 

} 

221 

222 

Ttol: :Ttol() { 

223 

FirstTol (0.1,0.1,100) ; 

224 

} 

225 

226 

void  Ttol: 

:FirstTol (double  prec,  double  resol,  double  Icmax) { 

227 

orgprecision  =  prec; 

228 

orgresolution  =  resol; 

229 

LCMAX=  Icmax; 

230 

InitTol ( ) ; 

231 

} 

232 

233 

void  Ttol: 

:  InitTol ( )  { 

234 

precision  =  orgprecision; 

235 

resolution  =  orgresolution; 

236 

} 

237 

238 

239 

Trigidice : 

:  Trigidice ( ) { 

240 

//  specific  gravity  of  ice? 

241 

// 

Y=rhoI/rhoW; 

242 

// 

lamda=l . Oe-6  ; 

243 

// 

eta=0 . 01 ; 

244 

// 

FW  =  1.0; 

245 

// 

FW  =  F(FW* (-rhoW*grav) ) ; 

246 

//  reduced 

latent  heat  of  fusion,  per  unit  volume  of  melt 

247 

// 

H  =rhoW*(lamda/gammaIA)*hIW; 

248 

InitTCALC 0 ; 

249 

} 

250 

251 

void  Trigidice :: InitTCALC 0 { 

252 

QW1=0; 

253 

QW2=0;//  reduced  mass  flux 

254 

DW=0;  //  change  in  water  content  across  layer 

255 

GRADUW1=0; 

256 

GRADUW2=0;//  reduced  water  pressure  gradient 

257 

QH1=0; 
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QH2=0;//  reduced  thermal  flux 
GRADTHETA1=0 ; 

GRADTHETA2=0 ; / /  reduced  thermal  gradient 
GRADPHI1=0; 

GRADPHI2-0 ; / /  reduced  gradient  of  PHI  from  Clapeyron 

Z=0;//  reduced  thickness 

ZN-0 ; //reduced  location  of  new  lens 

ZI=0 ; //reduced  location  of  old  lens 

DZ=0;//  reduced  layer  thickness 

UW=0 ;/ /reduced  water  pressure 

DUW=0;//  reduced  change  in  water  pressure 

//  across  layer 

UI=PHI;//  reduced  ice  pressure 

CHI=0 ;/ /reduced  value  of  the  stress  partition  factor 
UN=0;//  reduced  neutral  stress 

UNold=0;//  previous  value  of  reduced  neutral  stress 
UNmax=-1000 ; / /  maximum  value  of  neutral  stress 
THETA=0 ; 

DTHETA=0; 

GRADPROB=0; 

LAYERPROB=0 ; 

InitStress ( ) ; 

InitTol ( ) ; 

} 


void  Trigidice: : LAYER (double  rphi,  double  rdphi) { 

DW  =  W ( DSAT ( rphi ) ) -W (DSAT ( rphi+rdphi ) ) ; 

QW2  =  QWl  +  (VB-Y* (VB+VI) ) *DW; 

GRADUW2  =  FW  -  QW2 /RKW (rphi+rdphi ) ; 

QH2  =  QHl  +  H*(DW*VB  +  {QW2  ~  QWl ) ) ; 
GRADTHETA2  =  -QH2 /RKH (W (DSAT (rphi+rdphi ))) ; 
GRADPHI2  =  (Y-D*  GRADUW2  -  Y*H*  GRADTHETA2  ; 
if  ( {GRADPHI1<0)  I  I  (GRADPHI2<  0) )  { 

GRADPROB=l; 

DZ=1; 

}  else 

DZ  =  rdphi/ ( sqrt {GRADPHI1*GRADPHI2 ) ) ; 

DUW  =  FW*DZ  -  DZ* (QWl /RKW (rphi)  + 

QW2/ (RKW (rphi+rdphi) ) ) /2 ; 

DTHETA  =  ( (Y-1) *DUW-rdphi) / (Y*H) ; 

//  reset  variables  for  beginning  of  next  layer 
QW1=QW2 ; 

GRADUW1=GRADUW2 ; 

QH1=QH2; 

GRADTHETA1=GRADTHETA2 ; 

GRADPHI1=GRADPHI2 ; 

} 


void  Trigidice :: PROFILE () { 

LAYER ( PHI, DPHI) ; 
PHI=PHI+DPHI; 

DPHI  =  resolution/DZ ; 
THETA=THETA+DTHETA ; 
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Z=Z+DZ; 

UW=UW+DUW; 

UI=PHI+UW; 

CHI=SNYDER(PHI,DSAT(PHI) , DSAT ( PHI+DPHI ) ~DSAT { PHI ) ) ; 
UN=CHI*UW+ (l-CHI) *UI; 

} 

void  Trigidice :: GLOBAL {) { 

InitTCALC ( ) ; 

QWl  =  Y*VI  +  ( (Y-l) * (WSAT-WF) ) *VB; 

GRADUWl  =  FW  -  QWl/KWSAT; 

QHl  =  -RKH{WSAT) *GRADTB  ; 

FSI  =  VI/ (VI+VB) ; 

GRADTHETAl  =  GRADTB; 

GRADPHIl  =  (Y-D*  GRADUWl  -  Y*H*  GRADTHETAl; 

DPHI  =  resolution*PHIB*sqrt (GRADPHIl) ; 

DZ  =  DPHI /GRADPHI  Ir¬ 
resolution  =  DZ*DPHI; 

UW=0; 

UI=PHI; 

CHI=SNYDER(PHI,DSAT(PHI) , DSAT { PHI+DPHI ) -DSAT ( PHI ) ) ; 
UN^CKI *UW+ { 1 -CHI ) *UI ; 

THETA=-PHI/ (Y*H) ; 

} 


int  Trigidice : :prec 0 { 


if  {(  fabs ( (sum/ (ZI-ZN) -WF) /WD)  <=  precision) 
1  I  (LAYERPROB)  |  |  (GRADPROB)  )  { 

return  0 ; 

} 


else{ 


WF=sum/ (ZI-ZN) ; 


iteration++ 


return  1; 

} 

} 


int  Trigidice :: check  0 { 

if  (counter  >  LCMAX) { 

LAYERPROB=l; 
errors ( ) ; 
iteration  =  -999; 
return  0; 

} 

if  (GRADPROB) { 

errors ( ) ; 
iteration  =  -888; 
return  0; 
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366 

} 

367 

if 

(UN  >  UNmax) { 

368 

UNmax=UN; 

369 

ZN=Z, 

370 

return  1; 

371 

372 

if 

((UN  <=  UNmax) &&  (UI>UNmax) ) { 

373 

sum=sum+DW*DZ ; 

374 

ISI 

H 

II 

ISl 

375 

return  1 ; 

376 

377 

378 

if 

(  (UI<=UNmax) 

&&  (counter>=2) ) { 

379 

return  0; 

380 

} 

381 

return  1; 

382 

} 

383 

384 

385 

void 

Trigidice : : errors ( ) { 

386 

WF=0 

387 

G 

H 

II 

O 

388 

ZN=0 

389 

ISl 

H 

II 

1-^ 

390 

< 

H 

II 

o 

391 

} 

392 

393 

394 

void 

Trigidice : : once ( ) { 

395 

WF=WD; 

396 

PHI 

=  PHIB; 

397 

iteration  =  1; 

398 

do{ 

399 

counter  =  0 ; 

400 

PHI  =  PHIB; 

401 

sum  =  0 ; 

402 

GLOBAL ( ) ; 

403 

do  { 

404 

counter++ 

405 

PROFILE ( ) 

406 

}  while  ( 

407 

}  while  (  prec ( ) ) ; 

408 

} 

409 

410 

411 

412 

int 

Trigidice 

:pressdif f ( ) { 

413 

414 

415 

if 

(  f abs (Ulorg-UI)  <=  P (precision) ) 

416 

{ 

417 

return  0 ; 

418 

} 

419 

else  { 
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VI=VI+VIorg; 
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} 


// 


if  (/*  (LAYERPROB)  |  |  *  /  (UKUIorg)  } 

{ 

VI  =  VI  -  2.0*VIorg; 
VIorg=VIorg/2 . 0 ; 

VI  =  VI  +  VIorg; 
LAYERPROB  =  0; 

} 

if  (VKO.O) 

{ 

VI  =  V ( 0 . 5 /mmpdTomps ) ; 

} 


return  1; 

} 

int  Trigidice ; : slope { ) { 
if (GRADPROB) { 

orgresolution=orgresolution/2 ; 

LCMAX=LCMAX*2; 

return  1 ; 

} 

else  { 

return  0 ; 

} 

} 


void  Trigidice :: quasi () { 


} 


VB  =  VBorg; 

GRADTB  =  GRADTBorg; 

VI  =  V ( 0 . 5 /mmpdTomps ) ; 

VIorg  =  V(10/mmpdTomps) ; 

resolution=orgresolution; 

do{ 

do{ 

once ( ) ; 

}while (slope ( ) ) 
}  while  (pressdif f ( ) ) ; 


double  Trigidice :: heave () { 


quasi ( ) ; 

return  (VI  *  gammaIA*lamda/ (eta*nu) *  mmpdTomps) 

} 


double  Trigidice : :pressure () { 
once ( ) ; 

return  ( (UI*gammaIA/ lamda) /kP2P) ; 
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} 

double  Trigidice : : iterout { ) { 
return ( iteration) ; 

) 

double  Trigidice: : tgradout {) { 

QH2  =  -RKH(WSAT) *GRADTB+H* (VB* (WSAT-WF) +  Y*VI  +  ((Y-1)*(WSAT-WF))*VB); 

return{-l*theta0/eta*QH2*(FSI/RKHI()+{l-FSI) /RKH(WF) ) ) ; 

} 

double  Trigidice: :heatin{ ) { 

return ( -RKH (WSAT) *GRADTB/QH(1) ) ; 

} 

double  Trigidice : :heatout {) { 

QH2  =  -RKH (WSAT) *GRADTB+H* (VB* (WSAT-WF) +  Y*VI  +  ((Y-1)*{WSAT-WF))*VB); 

return {QH2/QH(1) ) ; 

} 

double  Trigidice :: water in () { 

QWl  =  Y*VI  +  {(Y-1)*(WSAT-WF))*VB; 

return (QWl*gammaIA*lamda/ {eta*nu) ) ; 

} 

double  Trigidice : : spacing ( ) { 

return ( (ZI-ZN) *VI/VB  *  eta  *1000); 

} 

double  Trigidice : : thickness ( ) { 

return(ZI  *  eta  *1000) ; 

} 


void  Trigidice :: look () { 

cout  << 

iteration  <<  " 

«  counter  «  "  ,  " 

<<  FSI  << 

<<  (UI*gaminaIA/lamda) /kP2P<< 

«  VI  *  gammaIA*lamda/ {eta*nu)  *  imnpdTomps  « 

«  VB  *  gainmaIA*lamda/ (eta*nu)  *  rampdTomps  « 

«  GRADTB*thetaO/eta  « 

QWl  =  Y*VI  +  { (Y-1) * (WSAT-WF) ) *VB; 

QHl  =  -RKH (WSAT) *GRADTB  ; 

QH2  =  QHl+H* (VB* (WSAT-WF) +QW1) ; 

//  fout  «  -1. 0*theta0/eta*QH2*  (VI*RKH(WF)  +VB*RKH(0)  )  /  (  (VI+VB)  *  (RKH  (  0 )  *RKH  (WF)  )  ) 

cout  «-1.0*theta0/eta*QH2*(FSI/RKHI()+(l-FSI) /RKH(WF) ) 

<<  "  ,  " 

<<  QWl*gainmaIA*lamda/ {eta*nu)  <<  " 

«  (ZI-ZN) *VI/VB  *  eta  *1000  « 

«  ZI  *  eta  *1000 
«  "\n"; 

} 
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HEAVE. CPP 


1  # include  "mcadincl .h" 

2  #include  "rigidice.h" 

3 

4 

5  LRESULT  RIGIDICEFunction (  LPCOMPLEXARRAY  calculated, 

6  LPCCOMPLEXARRAY  settings); 

7 


9 

10 

11 

12 

13 

14 

15 

16 

17 
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28 

29 
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FUNCTIONINFO  RIGIDICE  = 

{ 

"RIGIDICE",  //  Name  by  which  mathcad  will  recognize  the  function 

"settings",  //  heaverate  will  be  called  as  heavepressure ( settings ) 

"pressure,  heat  and  mass  fluxes  and  lens  sizes  for  given  settings",  // 

description  of  heavepressurre ( settings ) 

(LPCFUNCTION) RIGIDICEFunction,  //  pointer  to  the  executible  code 

COMPLEX_ARRAY ,  //  the  return  type  is  a  complex  scalar 

1,  //  the  function  takes  on  1  argument 

{  COMPLEX_ARRAY}  //  that  is  an  array 

}; 


LRESULT  RIGIDICEFunction ( 

{ 


LPCOMPLEXARRAY  calculated, 

LPCCOMPLEXARRAY  settings) 


Trigidice  foo; 
double  transient; 
double  micro,  macro,  body; 

double  khw,  khi,  khg; 

double  wsat,  wd,  ksat; 

double  alpha,  beta,  phib; 

double  prec ,  resol,  Icmax; 

double  heave,  pene,  gradtb; 


micro  =  settingS”>hReal [0 ] [ 0 ] ; 
macro  =  settings->hReal [0 ] [ 1] ; 
body  =  settings->hReal [0] [2]; 
f oo . FirstScales {  micro,  macro,  body); 


khw  =  settings -->hReal  [  0  ][  3  ]  ; 
khi  =  settings~>hReal [0]  [4]  ; 
khg  =  settings->hReal[0][5]; 
foo . InitTherm(  khw,  khi,  khg) ; 

wsat  =  settings->hReal [0] [6] ; 
wd  =  settings->hReal [0 ] [7 ] ; 
ksat  =  settings->hReal[0][8]; 
f oo . InitWATER {  wsat,  wd,  ksat) ; 

alpha  =  settings->hReal [ 0 ]  [ 9 ]  ; 
beta  =  settings->hReal [ 0]  [ 10 ]  ; 
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phib  =  settings->hReal [ 0] [11] ; 
foo.InitBC(  alpha,  beta,  phib) ; 

heave  =  settings->hReal [0]  [ 12 ]  ; 

pene  =  settings->hReal [0]  [13]  ; 

gradtb  =  settings->hReal[0][14]; 

foo . InitContinuum(  heave,  pene,  gradtb); 

prec  =  settings->hReal [0] [15] ; 
resol  =  settings->hReal [0] [16] ; 

Icmax  =  settings->hReal[0][17]; 
foo.FirstTol {  prec,  resol,  Icmax); 

MathcadArrayAllocate (  calculated,  1,  8,  TRUE,  FALSE); 

transient  =  foo .pressure () ; 
calculated->hReal  [0]  [0]  =  transient; 

transient  =  foo . tgradout ( ) ; 
calculated- >hReal  [0]  [1]  =  transient; 

transient  =  f oo . heatin ( ) ; 
calculated->hReal  [0]  [2]  =  transient; 

transient  =  foo . heatout ( ) ; 
calculated->hReal  [0]  [3]  =  transient; 

transient  =  f oo . waterin ( ) ; 
calculated->hReal  [0]  [4]  =  transient; 

transient  =  f oo . spacing {) ; 
calculated->hReal  [0]  [5]  =  transient; 

transient  =  f oo . thickness () ; 
calculated->hReal  [0]  [6]  =  transient; 

transient  =  foo . iterout { ) ; 
calculated->hReal  [0]  [7]  =  transient; 


return  0; 


//  return  0  to  indicate  there  was  no  error 


} 


LRESULT  heavepressureFunction (  LPCOMPLEXSCALAR  pressure, 

LPCCOMPLEXARRAY  settings); 


FUNCTIONINFO 

{ 

"heavepr assure" 
"settings" , 


heavepressure  = 


//  Name  by  which  mathcad  will  recognize  the  function 
//  heaverate  will  be  called  as  heavepressure ( settings ) 


"heave  pressure  for  given  settings" 

( LPCFUNCTION ) heavepressureFunction , 
COMPLEX_SCALAR , 

1, 

{  COMPLEX_ARRAY}  // 

}; 


,  //  description  of  heavepressurre (settings) 

//  pointer  to  the  executible  code 
//  the  return  type  is  a  complex  scalar 
//  the  function  takes  on  1  argument 
that  is  an  array 
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LRESULT  heavepressureFunction ( 


LPCOMPLEXSCALAR  pressure, 
LPCCOMPLEXARRAY  settings) 


Trigidice  foo; 
double  HeavePressure; 
double  micro,  macro,  body; 

double  khw,  khi,  khg; 

double  wsat,  wd,  ksat; 

double  alpha,  beta,  phib; 
double  prec,  resol,  Icmax; 

double  heave,  pene,  gradtb; 


micro  =  settings->hReal [ 0 ]  [0 ]  ; 
macro  =  settings->hReal [ 0 ]  [ 1 ]  ; 
body  =  settings->hReal[0][2]; 
f oo . FirstScales (  micro,  macro,  body); 

khw  =  settings->hReal [0]  [3 ]  ; 
khi  =  settings->hReal [0]  [4]  ; 
khg  =  settings->hReal[0][5]; 
f oo . InitTherm (  khw,  khi,  khg) ; 

wsat  =  settings->hReal [ 0 ] [ 6] ; 
wd  =  settings->hReal [0] [7] ; 
ksat  =  settings~>hReal [ 0 ]  [ 8 ]  ; 
foo . InitWATER(  wsat,  wd,  ksat)  ; 

alpha  =  settings->hReal [0]  [9]  ; 
beta  =  settings->hReal [ 0 ]  [ 10 ]  ; 
phib  =  settings->hReal [0 ] [ 11] ; 
foo.InitBC(  alpha,  beta,  phib) ; 

heave  =  settings->hReal [ 0 ] [ 12 ] ; 

pene  =  settings->hReal [0 ] [ 13 ] ; 

gradtb  =  settings->hReal[0][14]; 

foo . InitContinuum (  heave,  pene,  gradtb) 

prec  =  settings->hReal [0]  [15]  ; 
resol  =  settings->hReal [0 ]  [ 16 ]  ; 

Icmax  =  settings->hReal[0][17]; 
foo . FirstTol {  prec,  resol,  Icmax); 

HeavePressure  =  f oo . pressure () ; 
pressure->real  =  HeavePressure; 


return  0; 


//  return  0  to  indicate  there  was  no  error 


LRESULT  heaverateFunction ( 


LPCOMPLEXSCALAR  heaverate , 

LPCCOMPLEXARRAY  settings); 
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FUNCTIONINFO  heaverate  = 

{ 

"heaverate",  //  Name  by  which  mathcad  will  recognize  the  function 

"settings",  //  heaverate  will  be  called  as  heavepressure (settings ) 

"heave  rate  for  given  settings",  //  description  of  heavepressurre ( settings ) 
(LPCFUNCTION)heavepressureFunction,  //  pointer  to  the  executible  code 

COMPLEX_SCALAR,  //  the  return  type  is  a  complex  scalar 

//  the  function  takes  on  1  argument 
{  COMPLEX_ARRAY}  //  that  is  an  array 

}; 


LRESULT  heaverateFunction ( 

{ 


LPCOMPLEXSCALAR  heaverate , 

LPCCOMPLEXARRAY  settings) 


Trigidice  foo; 
double  HeaveRate; 
double  micro,  macro,  body; 

double  khw,  khi,  khg; 

double  wsat,  wd,  ksat; 

double  alpha,  beta,  phib; 

double  prec,  resol,  Icmax; 

double  press,  pene,  gradtb; 


micro  =  settings->hReal [0 ] [0] ; 
macro  -  settings->hReal [ 0 ] [ 1] ; 
body  =  settings->hReal [0] [2 ] ; 
foo.FirstScales (  micro,  macro,  body); 


khw  =  settingS“>hReal [0] [3 ] ; 
khi  =  settings->hReal [0] [4] ; 
khg  =  settings->hReal [0] [5]; 
f oo . InitTherm (  khw,  khi,  khg) ; 


wsat  =  settings->hReal [0] [6] ; 
wd  =  settings->hReal [0 ] [7 ] ; 
ksat  =  settings->hReal[0][8]; 
foo . InitWATER (  wsat,  wd,  ksat) ; 


alpha  =  settings“>hReal [0]  [9]  ; 
beta  =  settings->hReal [ 0]  [10]  ; 
phib  =  settings->hReal [0]  [11]  ; 
foo.InitBC{  alpha,  beta,  phib) ; 


press  =  settings->hReal [0] [12] ; 
pene  =  settings->hReal [ 0]  [13  ]  ; 
gradtb  =  settingS“>hReal [ 0] [14] ; 
foo , InitQuasi (  press,  pene,  gradtb); 


prec  =  settings->hReal [ 0] [15 ] ; 
resol  =  settingS“>hReal [ 0 ]  [ 16 ]  ; 
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Icmax  =  settings~>hReal[0][17]; 
foo . FirstTol (  prec,  resol,  Icmax); 

HeaveRate  =  foo . heave ( ) ; 
heaverate->real  =  HeaveRate; 


return  0; 


//  return  0  to  indicate  there  was  no  error 


} 

BOOL  WINAPI  DllEntryPoint  (HANDLE  hDLL,  DWORD  dwReason,  LPVOID  IpReserved) 

{ 

switch  (dwReason) 

{ 

case  DLL_PROCESS_ATTACH: 

{ 

//  DLL  is  attaching  to  the  address  space  of  the  current  process. 

// 

if  (  CreateUserFunction  (  hDLL,  ScRIGIDICE  )  =-  NULL  ) 
break; 

if  (  CreateUserFunction  {  hDLL,  Scheaverate  )  ==  NULL  ) 
break; 

if  (  CreateUserFunction {  hDLL,  &heavepressure  )  ==  NULL  ) 
break; 

CreateUserFunction (  hDLL,  &heavepressure  ); 


} 

// 


case  DLL_THREAD_ATTACH : 
case  DLL„THREAD_DETACH : 
case  DLL_PROCESS_DETACH: 

space . 


//A  new  thread  is  being  created  in  the  current  process. 
//A  thread  is  exiting  cleanly. 

//  The  calling  process  is  detaching  the  DLL  from  its  address 


247 

248  break; 

249  } 

250  return  TRUE; 

251  } 

252 

253 
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